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ABSTRACT 

We investigate the dynamics of putative Earth-mass planets in the habitable zone (HZ) of the 
extrasolar planetary system OGLE-2006-BLG-109L, a close analog of the Solar system. Our 
work is inspired by work of Malhotra and Minton (2008). Using the linear Laplace-Lagrange 
theory, they identified a strong secular resonance that may excite large eccentricity of orbits 
in the HZ. However, due to uncertain or unconstrained orbital parameters, the sub-system of 
Jupiters may be found in dynamically active region of the phase space spanned by low-order 
mean-motion resonances. To generalize this secular model, we construct a semi-analytical 
averaging method in terms of the restricted problem. The secular orbits of large planets are 
approximated by numerically averaged osculating elements. They are used to calculate the 
mean orbits of terrestrial planets by means of a high-order analytic secular theory developed 
in our previous works. We found regions in the parameter space of the problem in which 
stable, quasi-circular orbits in the HZ are permitted. The excitation of eccentricity in the HZ 
strongly depends on the apsidal angle of jovian orbits. For some combinations of that angle, 
eccentricities and semi-major axes consistent with the observations, a terrestrial planet may 
survive in low eccentric orbits. We also study the effect of post-Newtonian gravity correction 
on the innermost secular resonance. 

Key words: celestial mechanics - secular dynamics - relativistic effects - analytical methods 
- extrasolar planetary systems 



1 INTRODUCTION 

Recently, |Gaudi et aL] ( |2008| ) announced a discovery of two-planet 
extrasolar system hosted by the OGLE-2006-BLG-109L star of 
~ 0.5 Mq. Its jovian companions of ~ 0.71 and ~ 0.27 Jupiter 
masses are in orbits of 2.3 and 4.6 au, resembling a scaled copy of 
the Jupiter- Saturn pair. This similarity leads to a natural question, 
whether additional planets can exist in this system, in particular 
Earth-like planets in the habitable zone (HZ). As for now, the an- 
swer can be speculative, because the OGLE-2006-BLG-109L plan- 
etary system has been detected during an observed micro-lensing 
event and hence the observational data are not reproducible. In ad- 
dition, the faint and distant parent star cannot yet be observed by 
other planet-hunting techniques, like the Doppler spectroscopy or 
astrometry, to eventually detect the presence of an Earth-like planet. 
Moreover, some orbital parameters of the jovian planets are either 
uncertain or undetermined at all. The orbital architecture of the sys- 
tem, including putative terrestrial planets, and its stability, may be 
investigated only indirectly by numerical simulations. 

In this work, we focus on the dynamical structure of the HZ. 
The HZ is defined implicitly through a requirement of the liquid 
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state of water (see, e.g., |Kasting et al.|1993||Hinse et al.|2008[|Mal-| 
[hotra & Minton|20'08] ). That condition is satisfied if the terrestrial 
orbit, parameterized by the semi-major axis ao and eccentricity eo, 
is confined to the annulus of {rp^Va) = (0.25,0.36) au [following 
[Hinse et ar|p008| ), the terrestrial planet in the HZ will be called 
with Latin word Tellus from hereafter and we will always refer to 
its orbital elements with index "0"]. It means that the pericenter and 
apocenter distances of habitable orbit must be limited through the 
following conditions : 

ao{\-eo)>rp and ao(l+^o)<^a, (1) 

providing very long time- scale of stable motion, well enough to 
develop and support life (see, e.g. |Hinse et ar]|2008| for details). 
The term stable or HZ-stable, regarding the terrestrial orbit, will 
be understood as "confined to the HZ, as defined above, during the 
secular time scale counted in tens of Myrs". Apparently, the HZ 
of the OGLE-2006-BLG-109L system is separated from the res- 
onant influence of the primaries. The ratio of orbital periods cor- 
responding to ao G (0.25,0.36) au is ~ 20. Therefore, no strong 
mean motion resonances (MMRs) can perturb putative terrestrial 
orbits. However, even in the absence of the MMRs, such orbits 
may be still disturbed by the long-term, secular interactions with 
jovian companions. Indeed, [Malhotra & Minton| ( |2008| ) identified 
two secular resonances (Vi^2) in the inner part of the OGLE-2006- 
BLG-109L system related to the fundamental secular frequencies 
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(modes) of the orbits of Jupiters. These resonances are analogs of 
the V5 6 resonances with Jupiter and Saturn in the Solar system. 
Both of them lead to strong amplification of ^0 to large values (up 
to 0.8-0.9), and are centered around ao ~ 0.3 au, and ao ~ 0.7 au, 
respectively. In particular, the Vi resonance acts in the region of 
HZ, suppressing long-term stable motion confined to that zone. The 
motion of a terrestrial planet may be stabilized through a change of 
the secular frequencies by interactions with an additional (yet un- 
known or hypothetical) inner planet ( [Malhotra & Mintonl|2008| ). 
Such an object of ~ 0.3 and semi-major axis ~ 0.1 au would 
change the Vi secular frequency, and shift the "troublesome" reso- 
nance out of the HZ. Similarly, numerical studies of |Pilat-Lohinger| 
|et al.| ( |2008| ); |Pilat-Lohinger et al.| ( ^2008j , regarding Solar-like sys- 
tems show that a small planet in the inner region may change the 
position of the secular resonance V5 and significantly affect the 
habitability for a given semi-major axis. Earlier, the long term ef- 
fects of adding or removing planets in the inner Solar system were 
analyzed by ,Namouni & Murray, ( (1999j ), following [Innanen et al.| 
( |1998| ) who found that the Earth-Moon (EM) stabiHzes the orbits 
of Venus and Mercury by suppressing a strong secular resonance 
of period 8.1 Myr near Venus orbit. Actually, Namouni & Murray] 
( |1999| ) conclude that is not clear whether the absence of EM does 
really provide a physical insight into the stability of the whole sys- 
tem. They argue that throughout the accretion ofplanetesimals, the 
mass that formed Earth contributed consistently to the secular pre- 
cession rates and more generally to the current state of the Solar 
system. Hence, the studies of long-term behavior of a planetary sys- 
tem with arbitrary initial conditions (e.g., modified by adding or re- 
moving planets while keeping other orbital elements fixed), which 
are partially consistent with observations, seem not well justified 
in general. Hence, because the OGLE-2006-BLG-109L- system ar- 
chitecture is not fully determined (in fact, it is known poorly, as 
we will see below), rather to analyse the influence of a few hy- 
pothetical bodies in the system on the dynamical structure of the 
HZ, we follow the known observational constraints on the system. 
Its Jovian planets are orbitally coupled, and the habitable zone is 
very close to the star, prohibiting extensive numerical integrations. 
Hence, the main motivation of our work is to build up and inves- 
tigate an appropriate semi- analytic, specific model of the system, 
and try to answer whether putative, single Earth-like planet may 
survive in the HZ during the secular time- scale. 

In this work, we are heavily inspired by the ideas and results 
of [Malhotra & Minton| ( 2008). To model the secular dynamics of 
the OGLE-2006-BLG-109L system, they appUed the classic lin- 
ear Laplace-Lagrange (L-L) theory (see, e.g., [Murray & Dermott] 
[200 0), assuming non-resonant, quasi-circular orbits of the jovian 
planets. Yet, as we shall demonstrate in Sect. 2, the orbital data 
combined with the best-fit errors permit a variety of orbital states, 
including resonant or highly eccentric configurations. In such cir- 
cumstances, the L-L theory cannot give proper estimates of the 
secular frequencies. In this work, we developed a relatively simple 
quasi- analytic algorithm that helps us to avoid the limitations of the 
L-L theory. We focus on the global dynamical structure of the HZ, 
taking into account wide ranges of the orbital parameters of jovian 
companions permitted by the best-fit uncertainties. In particular, we 
analyse the sensitive dependence of the secular motion in the HZ on 
unconstrained apsidal angle AG5b,c = G5b — G5c, where G5b and G5c are 
the longitudes of periastra of planets b and c, respectively. We also 
investigate the long-term stability of the OGLE-2006-BLG-109L 
system by means of dynamical maps constructed with the fast indi- 
cator MEGNO, (7), (Cincotta & Sim6|2000|[Cincotta et al.|2003 
and direct numerical integrations. 



2 DYNAMICS OF THE OGLE-2006-BLG-109L SYSTEM 



According to the discovery paper ([Gaud i et al.|2008| ), the orbits of 
jovian companions of OGLE-2006-BLG-109L are determined with 
significant uncertainties, i.e., = 2.3 ±0.2 au, a^ =4. 6 ±0.5 aufor 
the inner and outer planet, respectively; ^b is unconstrained at all, 
and Cc = 0. 1 1 ± 0.07. To keep our work consistent with a parame- 
terization of the analytical theory in ( Migaszewski & Gozdziewski] 
[2008| ), we interpret these orbital parameters as the canonical, ge- 
ometric elements related to the Poincare coordinates (see, e.g., 
[Morbidelli[[2002l [Ferraz-Mello et al.[[2006l ). From the "practical" 
point of view, these elements are not significantly different from 
the common, astrocentric Keplerian elements. The longitudes of 
pericenters (G5b,G5c), and mean anomalies are undetermined. The 
inclination of OGLE-2006-BLG-109Lc has been estimated from 
observations as / ~ 59°, and also masses of companions are con- 
strained within ~ 10% error. It is reasonable to assume that the 
whole system is coplanar. Still, significant errors of the semi-major 
axes leave room for a few qualitatively different orbital configura- 
tions which can be studied with the help of 2D dynamical maps. We 
select and vary two Keplerian elements, and other orbital parame- 
ters are fixed. For each point of the parameter plane, representing 
an initial configuration of the system, we compute the MEGNO 
signature ( Cincotta & Sim6[[2000| ) with the symplectic algorithm 
( [Gozdz iewski et al.[[2008| ). A relatively short integration time re- 
quired by MEGNO and good sensitivity of the indicator to chaotic 
motions make it possible to investigate large volumes of the phase- 
space with high resolution, and to detect efficiently unstable con- 
figurations. For instance, the dynamical maps in the (ab,^b) -plane, 
illustrated in the left-hand panels of Fig. [T] have the resolution of 
640 X 400 points and each point in these map represents a config- 
uration integrated over ~ 0.4 Myr. The orbital parameters of the 
nominal configuration marked with crossed circle in the map are 
taken from (Mal hotra & Minton|2008| ), and we choose AG5b,c, fol- 
lowing a concept of the so called representative plane of initial 
conditions (Michtchenko & Malhotra 2004). That (^b cos AG5b,c, ^c)- 
plane, S from hereafter, is defined through fixed AG5b,c = (the 
right-hand half-plane, for initially aligned orbits) and for AG5b,c = 7i 
(orbits are anti- aligned). For these particular values of AG5b,c, the 
eccentricities in two-planet configuration reach maximal or mini- 
mal values at the same time with an accord to the conservation of 
the total angular momentum and the condition of ^i^ec/^ AG5b,c = 
is for the secular Hamiltonian of the three-body problem). It 
can be shown ( Michtchenko & Malhotra 2004) that all phase tra- 
jectories of non-resonant two-planet system must intersect the 5- 
plane. Hence, the global dynamics of the system may be described 
conveniently in terms of initial conditions selected in the 5-plane. 
The 5-plane is particularly useful for investigating equilibria of the 
secular system. For a reference, other examples of MEGNO maps 
in the 5-plane are shown in the right-hand panels of Figs.[T]and[2] 

Clearly, the pair of Jupiters resides in dynamically active re- 
gion of the phase space which is spanned by a few low-order 
MMRs. Allowing for a variation of ab within the 0.2 au error, 
the OGLE-2006-BLG-109L system can be found in the neighbor- 
hood of the 3c: lb MMR or 5c:2b MMR. Because the error of a^ is 
~ 0.5 au, even 2c: lb MMR configurations are permitted. Simulta- 
neously, as the dynamical maps indicate, the border of stable zone 
strongly depends on initial eccentricities. Therefore, to determine 
the secular evolution of orbits in the HZ, stable configurations of 
the giant planets must be chosen with care. To illustrate that is- 
sue, we select a few representative configurations (still, consistent 
with the observational uncertainties) and compute their dynamical 
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Figure 1. Dynamical maps of the OGLE-2006-BLG-109L system in terms of the MEGNO indicator, (7). The type of orbits is color-coded; yellow means 
strongly unstable configurations of jovian companions, black means quasi-periodic solutions. The integration time is ~ 10^ orbital periods of planet c (~ 
0.4 Myr). The position of the nominal system in the phase space is marked by crossed circle in each map. Left-hand panels are for the (<2b,^b) -plane and 
A03b,c = 0,71, respectively. Panels in the right-hand column are for the representative plane (see the text for details) and the MEGNO range set to [0,100] (the 
top panel) and [1.98,2.18] (the bottom panel). A comparison of these maps reveals some dynamical structures and transition zones between strongly and mildly 
chaotic motions. Note a change of the color coding in the bottom right-hand panel. 



Table 1. Orbital models of the OGLE-2006-BLG-109L system considered 
in this work. Masses of Jovian planets are fixed to 0.71 mj and 0.27 mj 
for planets b and c, respectively, following Gaudi et al. (2008). Mass of 
the parent star is 0.5 mq. The mean anomalies fA4,c in these models are 
constrained by the mean longitude A,b,c = ACj3b,c + ^h,c = 0- 



Model 


«b [au] 


flc [au] 


eh 


ec 


Acnb,c 


k 


2.30 


4.60 


0.03 


0.05 





In 


2.30 


4.60 


0.03 


0.05 


71 


Ik 


2.30 


4.60 


0.10 


0.10 





IIIo 


2.40 


4.35 


0.10 


0.10 





Illn 


2.40 


4.35 


0.10 


0.10 


71 


IVo 


2.10 


5.10 


0.10 


0.10 





IVn 


2.10 


5.10 


0.10 


0.10 


71 


Vo 


2.50 


4.10 


0.45 


0.10 






maps. Orbital parameters of these models, in terms of parameter 
tuples (ab [au] , [au] , ^b , , AG5b,c ) . are given in Table [l] model /p 
for the nominal configuration investigated in ( [IMalhotra & JVIinton] 
|2008| ), with apsides aligned and model /j^ with apsides anti-aligned; 
model //q with moderate eccentricities, close to the best fit elements 
quoted in (Ga udi et al.,,20Q8 ); models ///q and ///t^, for compact 
systems with moderate eccentricities; models /V(3,7i for hierarchical 
configuration of the Jupiter s, and model Vq involving Jupiter s in the 
2c: lb MMR. Dynamical maps for models ///o,7r are shown in the 
left-hand panel of Fig. [2] and for model Vq in two remaining panels. 



3 SEMI- ANALYTIC MODEL OF THE HZ 

The time-scales of orbital periods of giants and Tellus in the HZ 
are very different. Moreover, due to very long period of the Vi res- 



onance, of a few Myrs (~ 10^ orbital periods of Tellus), the di- 
rect integrations of the planetary equations of motion would re- 
quire incredible amounts of CPU time. Therefore, we consider Tel- 
lus placed in the HZ and each planet in the dynamically coupled 
pair of Jupiters as highly hierarchical system, because the ratio of 
semi-major axes a^ja^^c ^ 1/10 is small. That makes it possible to 
average out the short-term variations of the orbit of Tellus over the 
mean longitude, and to approximate its secular evolution by means 
of a 24-order analytical theory in the semi-major axes ratio ( |Mi-| 
[gaszewski & Gozdziewski|2008| ). Yet the theory cannot be directly 
applied to jovian orbits, because their elements can be varied in 
wide ranges, and non-resonant as well as resonant or near-resonant 
configurations are permitted. The averaging of the whole system, 
with strongly interacting giant planets cannot be done uniformly 
(see, for instance [Malhotra et al.|1989| l. Hence, we introduce two 
simplifications to the problem. At first, we assume that Tellus does 
not affect the orbital evolution of the giants, so we consider the re- 
stricted problem. Then the secular Hamiltonian (per mass unit) may 
be written as a sum of two terms i^ec = + , for interactions 
of Tellus with each Jupiter, separately. 




where ao i = a{)/ 

AO; 



oco,/ 



^^^''■^(^o,^/,AG5o,-) 



and / = 



b,c. The explicit formulae for expan 

r 



sion terms, (g o,g/,AG3oj-) can be found in (Migaszews ki &| 

Gozdziewski|2008| ) . The equations of motion implied by i^ec re- 



fer to the mean orbital elements of Tellus and the mean orbital 
elements of Jupiters as parameters of the model. In the realm of 
the restricted problem, we can derive the mean orbits of Jupiters 
(and their perturbations) by numerical averaging. The full, N-bo&y 
equations of motion of these planets are integrated numerically at 
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Figure 2. Dynamical MEGNO maps for two putative configurations of the OGLE-2006-BLG-109L system. The left-hand, and the middle panels are for the 
5-plane (see the text for details), the right-hand panel is for the (flbj^b) -plane. Orbital elements of tested configurations are marked with crossed circle and 
labeled. The left hand panel is for model ///, the middle- and the right-hand panels are for the 2c: lb MMR of Jupiters (model Vq)- 
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time [ky] time [ky] semi-major axis [au] 

Figure 3. The semi-analytic theory of Tellus explained in graphic form. Panel a). A part of ec{t)-cm\Q (dots, the full three-body numerical integration), and 
the averaged eccentricity of the outer planet (red, filled circles). Panel b). The mean eccentricity of the outer planet after averaging. Panel c). A comparison of 
models of the Vi secular resonance derived with the help of the semi- analytical theory (smooth, black curve) and with the direct numerical integration of the 
restricted four-body problem (open circles). Orbital parameters of the jovian sub-system are written in the box (model IV-n). Orbits with largest max^o ^ 0.8-0.9 
are chaotic and that explains significant deviations of the numerical values from the semi-analytic solution. The red double curve marks geometric boundary 
of the HZ impHed by Eq.[T] Initial conditions are: A03o,c = 0, = 0.01, and the mean anomahes of Tellus were chosen at random for the numerical 
integrations. 



least over of a few Vi periods (~ 10 Myr). Then the temporal, 
osculating elements are converted with appropriate constant time 
step- size to the mean elements through the running average (see 
Figs.[3^,b). Because we intent to integrate the secular equations of 
motion of Tellus, the time step of the averaging (- 10^ yr) should 
be a small part of the apsidal period of its orbit. With this relatively 
long time step, as compared to the orbital periods, we remove all 
fast, quasi-periodic variations of the osculating elements. The mean 
orbits of jovian planets need to be computed only once (because 
we consider the restricted problem), and then we can efficiently 
reconstruct secular orbits of arbitrary number of Tellus "clones", 
integrating numerically the secular equations of motion induced by 
i^ec- The integration of a single terrestrial orbit over the secular 
time scale (typically, a few Myrs) is rapid, by a few orders of mag- 
nitude shorter than of the full system: the CPU time of the analytical 
solutions counts in minutes, but to integrate ~ 100 orbits of Tellus 
during - 20 Myrs, the MERCURY code ( Chamber s|p^99l inte- 
grator RADAU) spent a week on four 2GHz AMD/Opteron CPU 
cores. A test of this approach is illustrated in FigJsL According 
with the classic L-L theory ( [Murray & Dermottl|2000| ), the cen- 
tral peak of eccentricity eo appears when the particle's apsidal fre- 
quency is close to the forcing planetary frequency. Hence that peak 
shows the position of the resonance in the parameter space (for 
instance, at the ao-axis). The position of the eccentricity peak as 
well as the shape of max^o(^^o) graph are reproduced with great 



precision in the whole range of ao, in spite of extremely large ec- 
centricity attained by terrestrial orbits around ao ~ 0.40 au. Actu- 
ally, after a few Myrs such orbits may become strongly chaotic. 
We did similar tests, choosing angle AG5o,c and the mean anomaly 
at random. The results are illustrated in Fig. [4] At this time, 
the orbits of Jupiters are set as in models I-V (see Table [TJ. We 
compute the maximal eo of Tellus for varied initial ao in the rele- 
vant range of [0.2,0.4] au. We selected ~ 100 of clones of Tellus, 
setting their initial eo = 0.0\. Next, we compared max^o attained 
during 10 Myrs, as computed with the help of the quasi- analytic 
algorithm and during 10-20 Myrs of the numerical integration of 
the restricted four-body problem. Figure |4] is for one-dimensional 
plots of max^o as a function of initial aQ. Smooth curves obtained 
for AG5o,c = G5o — G5c = 0,71 are for the semi-analytic theory, filled 
circles are for the numerical integrations. Elements of Jupiters are 
labeled in boxes drawn in each respective panel and with an appro- 
priate Roman number referring to Table[T] The results for all tested 
configurations perfectly coincide. Filled points are found strictly in 
the limits by analytical curves obtained for initial orbits of Tellus 
aligned, and anti- aligned with planet b (or planet c). The position 
of the Vi resonance, the width of the eccentricity peak, and the 
maximal range of eo are determined with great precision. These 
data are critical for estimating a fraction of stable, quasi-circular 
orbits in the HZ: for a given initial ao we compute the secular 
value of max^o and then we compare that value with the limits 
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Figure 4. The maximal eccentricity of Tellus as a function of initial semi-major axis ao. Initial eccentricity eo = 0.01. Mass parameters of the OGLE-2006- 
BLG-109L system are: mo = 0.5 M©, mb = 0.71 mj, mc = 0.27 mj. Subsequent panels are for initial orbital elements of the Jupiters expressed in terms of 
parameter tuples: (^b [au] , ac [au] , ^b , , A03b,c ) labeled in the boxes (see also the text and Table[T]for details). Smooth grey /black curves are derived through the 
semi-analytic theory (for initial A03o,c = 0, tt, respectively), filled circles are for the results of the direct numerical integrations spanning ~ 20 Myrs for random 
A03o,c- Green areas mark ranges of ao providing stable orbits of Tellus confined to the annulus of HZ between 0.25 au and 0.36 au for A03o,c = 0- 



(Eq. [TJ implied by the orbital annulus of the HZ. Its geometrical 
borders are marked with red solid curves in Fig. |4] also the rele- 
vant range of semi-major axes, providing secularly low-eccentric 
HZ-orbits for initial AG3o,c = 0» are marked by projections of green 
areas under the max^o curve onto the ao-axis. Panels of Fig.|4^,b 
are for the configuration studied by |Malhotra & M inton ( 2008 ), i.e., 
e\) = 0.03, = 0.05, and ab = 2.3 au, ac = 4.6 au. Panel of Fig.|4^ 
is for aligned orbits (our model /q), while Fig.|4]3 is for anti-aligned 
orbits (our model Iji). We can observe significant differences be- 
tween these plots, in spite of that initial eccentricities and semi- 
major axes are the same, and only initial relative orbital phase of 
jovian planets (in fact, AG5b,c has changed). That is quite a different 
conclusion than could be derived in the framework of the L-L the- 
ory. For AG5b,c = a large part of the HZ is rendered unstable by 
the secular resonance, nevertheless, for the anti-aligned configura- 
tion, most particles would survive in the HZ, and max^o(^o) would 
be relatively small. Moreover, the instability generated by the sec- 
ular resonance (with the centre at approximately 0.29 au) is much 
stronger when initial eccentricities of Jupiters are ~ 0.1 (Fig.|4]:). 
In such a case, max^o can be as large as 0.8-0.9 and all telluric 
orbits will be wandering far out of the HZ. Similarly, for moderate 
^b,c but different semi-major axes of the jovian system (Fig.|4|i,e,f), 
the position of the secular resonance can be significantly shifted. 
Still, there will be some ranges of ao permitting secularly stable 
motion of Tellus in the HZ (we recall that these ranges are marked 
in green). A particularly interesting and quite surprising result is 
illustrated in the last panel of Fig.|4f, that refers to the 2c: lb MMR 
of the jovian planets. This MMR would protect Tellus from the sec- 
ular excitation of the eccentricity, through shifts of the max eo peak 
generated by Vi far out of the HZ. Also this case confirms excellent 
accuracy of the semi-analytic approach. As we can observe, even 
if the Vi resonance disturbs the orbits of Tellus, usually, there is a 
range of ao, Aao, for which its orbit will remain entirely in the HZ. 



According with the definition of the HZ, the maximum range of ao 
is maxAao ~ 0.1 1 au. To measure a fraction of HZ-stable orbits as 
a function of ao, with respect to different configuration of jovian 
planets, we define the linear coefficient of habitability: 



/h: 



Aao 
maxAao ' 



/hzG[0,1] 



(2) 



where Aao is the range of initial semi-major axis implying orbits 
entirely confined to the HZ during the secular time-scale of a few 
Myrs. If the whole HZ is rendered unstable by the Vi resonance 
(or other perturbation) then fuz = 0, and if all orbits remain within 
the HZ annulus, we have fuz = 1- For instance, for the nominal 
system in Fig. [4^, /hz ~ 0.4. The next configuration illustrated in 
Fig.|4|D has fuz ~ 0.95; in this case, the very narrow Vi resonance 
affects a small part of the HZ. We can observe in Fig. |4] that the 
position and the width of Vi depends on assumed orbital parameters 
of the jovian sub- system, hence fuz = /Hz(^2b,<3c,^b,^c5 AG5b,c), also 
parameterized by masses of large bodies mo, mb, mc. In general, fuz 
depends also on initial eo and the relative phase of Tellus and a 
jovian planet, e.g., AG5o,c. 

We also investigated the influence of the general relativity 
(GR) correction on the secular dynamics of Tellus. In our recent 
paper (Migaszewski & Gozdziewski 2009 ) we have shown that 
this apparently subtle correction to the Newtonian gravity may af- 
fect the secular dynamics significantly, and they can be particularly 
important for small and close-in planetary companions. With the 
semi-analytic model of the jovian orbits, we could repeat the cal- 
culations of max^o^ adding the GR correction to i^ec- The results 
derived with these corrections (the relativistic model) and with the 
Newtonian interactions only {the classic model) are compared in 
Fig.|5] For the nominal configuration, we observe only a small shift 
of the eccentricity peak (of ~ 0.01 au). The maximal eccentric- 
ity is not affected significantly. Indeed, the GR corrections become 
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important when the rate of the apsidal motion, which it causes to 
change, is similar to the effect of the Newtonian point-mass interac- 
tions. In the region of HZ, the GR induced apsidal advance has the 
period of a few Myrs, while the Newtonian period is ~ 30,000 yr 
only. Hence, we may conclude that the GR effects are not really 
important for the stable motion of Tellus in the HZ, and we do not 
consider them anymore. 



4 SIZE OF THE HZ IN THE PARAMETRIC SPACE 

The one-dimensional parametric survey of the HZ with respect to 
aQ (Fig.[4| can be generalized to the second dimension of the initial 
eccentricity cq. Similarly to /hz, we define the planar coefficient of 
habitability: 

^0 

■^HZ ='^Hz(^«b,^2c,^b,^c,AG5b,c) = ^HZ ^ [0,1], (3) 

maxEo 

where maxEo is the area of initial conditions in the given parameter 
plane consistent with the geometrical boundaries of the HZ, while 
Eq measures the area of initial conditions in that plane, implying 
orbits confined to the HZ during the secular time-scale. In contrary 
to /hz, the planar coefficient of habitability does not depend on cq. 
Fixing AG5o,c = 0, 7i, we follow again the concept of the representa- 
tive plane of initial conditions. 

The results for the same configurations of Jupiters as in Fig.|4] 
are illustrated in Fig. [6] Because the stable zone is bordered by or- 
bits of Tellus initially aligned or anti- aligned with the apsidal line 
of planet c (or, alternatively, with planet b), both these limiting 
cases can be conveniently shown in one (ao,^o cos AG5o,c) -plane. 
The negative values on the j-axis tell us that AG5o,c = 7i, positive 
values mean that AG5o,c = 0. The borders of the zone with HZ- 
stable orbits are marked with green solid curves. The geometric 
boundary of maxEo is marked with thick dotted curves. The stable 
zone has rather complex shape, and, depending on orbits of the jo- 
vian sub-system, we can derive quite unexpected conclusions. For 
their apsides aligned (Fig.|6^), there is still significant area of hab- 
itable orbits for initial cq ^ 0. Curiously, if the jovian orbits are 
initially anti-aligned, the HZ is almost wholly preserved (Fig.[6|D). 
In the next instance (panel c), the HZ is rendered unstable, while 
in all remaining cases (Fig.[6|i,e,f), there are HZ-stable orbits up to 
eo ~ 0.15. 

Yet the tests illustrated in Figs. |4] and [6] still provide limited 
information on the HZ because we analyzed only isolated configu- 
rations of Jupiters. To obtain better insight into the dependence of 
the HZ on different orbital setups of Jupiters, we extend the survey 
to computation of ^hz in two-dimensional parameter planes of the 
jovian sub-system. 



4.1 The HZ in the representative plane 

The features of the 5-plane make it particularly convenient to anal- 
yse suz for fixed semi-major axes of the jovian sub-system and its 
dependence on the eccentricities. The results are illustrated in the 
left-hand panel of Fig. [7] which shows color-coded values of suz 
in the 5-plane constructed for the nominal system (ab = 23 au, 
ac = 4.6 au, our model /). It can be regarded as non-resonant or 
near-resonant (see Fig. [T]). In such a case, in the regime of mod- 
erate eccentricities, the apsidal angle AG5b,c circulates or librates 
around or 7i (these regions are dotted in Fig. [7]). For a refer- 
ence, the red, thick curves mark approximate positions of stable 



stationary solutions of the non-resonant secular model of the jo- 
vian sub-system, calculated with the help of a high-order analytic 
theory (Migaszewski & Gozdziewski 2008 ) . The curve in the pos- 
itive half-plane of S is for the so called mode I equilibria, charac- 
terized by librations of AG5b,c around for close phase trajectories; 
mode II curve (in the negative half-plane of S) is surrounded by or- 
bits with secular angle AG5b,c librating around n [see ( Michtch enkoj 
|& Malhotra|2004t ) for details]. The right-hand panel of Fig.|7]is for 
the phase diagram constructed for a fixed angular momentum in 
the (^bC0sAG5b,c,^b sin AG5b,c) -plane, illustrating the equilibria and 
neighboring phase trajectories, including the nominal configuration 
marked with open circles. White regions in the left-hand panel of 
Fig.[7]indicate eccentricities and AG5b,c for which the HZ is strongly 
affected by the Vi resonance, hence suz = 0. Black regions are for 
suz ~ 1- Before integrating secular orbits of Tellus, we must av- 
erage out orbits of primaries, hence we have also "an occasion" 
to eliminate unstable configurations (disrupted during ~ 10 Myrs, 
(Y) > 10). They are marked with gray crosses. Nominal eccentric- 
ities of the Jupiters are marked with filled circles (labeled with Iq 
for AG5b,c = 0, and Ijz for AG5b,c = For the right-hand half-plane, 
Suz ~ 0.5, and for the left-hand half-plane of 5, suz ~ 0.95. Of 
course, it must be in accord with the results of two-dimensional 
surveys (see Fig.[6| performed for the same initial conditions. 

Clearly, the region of large ^hz lies in the neighborhood of 
stationary mode II. On contrary, the secular resonance in the con- 
figurations selected in the positive half-plane (AG5b,c = 0) preserves 
stable orbits in a very small part of the HZ. This seems a general 
feature of the system, and in fact it can be explained through the 
fundamental frequencies. Close to the equilibria, AG5b,c oscillates 
around 0, or n. To the first significant octupole terms of i^ec, 

deo 

— — ~ AsinAG5ob +5sinAG5oc, 
at 

where A and B are coefficients dependent on the orbital elements, 
having the same sign. Hence, the anti-aligned orbits of the giant 
planets favor small (or much smaller) excitation of cq than the 
aligned jovian orbits. 

The results of the next simulation are shown in Fig. [8] Now, 
we changed the semi-major axes to ab = 2.5 au, ac =4.1 au (our 
model Vq, in the 2c: lb MMR). We recall that the middle- and the 
right-hand panels of Fig. |2] are for the dynamical maps of this 
configuration in the (^b,^c)-, and (ab,^b) -planes. Figure [s] shows 
the results in terms of color-coded Suz plotted in the right-hand 
half-plane of S (there is no equivalent of a stable region in the 
left-hand half-plane of S). We compute the suz for jovian orbits 
stable at least over ~ 10 Myrs. There is a perfect match of the 
border of stable zone derived through the long-term integrations, 
and with the MEGNO map (both maps are plotted in Fig. [8] see 
also Fig. [2]), although the indicator was calculated over ~ 0.4 Myr 
"only". Comparing these two maps, we found a transient zone 
around (abC0sAG3b,^b) ~ (0.45 au,0.05), in which MEGNO al- 
ready detects chaotic motions, while the system is disrupted after a 
few Myrs. It confirms a precise calibration of the integration time 
of(F>. 

4.2 The habitable zone in the (ab,flc) -plane 

Finally, we computed ^hz in the (ab,ac) -plane, fixing the reference 
values of ^b = 0.03, = 0.05 (our model /). For a comparison, we 
choose two values of AG5b,c = 0,7i. The right-hand panel of Fig.|9] 
is for initially aligned orbits of Jupiters, the left-hand panel is for 
initially anti-aligned orbits. The position of the nominal system in 
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Figure 5. A comparison of the maximal eccentricity of Tellus derived with the help of the semi-analytical theory, in the range of relevant for the HZ in 
model /. The initial eccentricity cq = 0.01 and A03o,c = 0. The left-hand panel is for Atnb,c = 0. the right-hand panel is for A03b,c = ^- The black curve is for 
the model with the general relativity corrections, the gray curve is for the classic, Newtonian point-mass model. 
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Figure 6. Color coded maximal eccentricity of Tellus in the (ao^^o cos ACj3o,c) -plane; A03o,c = in the positive half-plane, and A03o,c = tt in the negative 
half-plane. The thick, dotted curves mark geometric borders of the HZ with accord to condition (Eq.[^ see the text for details). The green solid curves border 
initial conditions of Tellus implying orbits confined to the HZ during at least ~ 10 Myrs. Subsequent panels are for the orbital configurations of Jupiters given 
in Table[T]and labeled accordingly, see also Fig. |4] (models I-V). 



the parameter space is marked with filled circles. Crosses, along 
straight lines, indicate unstable solutions detected with the help of 
MEGNO and, in fact, they are related to the relevant, low-order 
MMRs, e.g., 5c:2b, 3c: lb and 7c:2b. Clearly, the jovian system 
would be more "friendly" for inhabitants of Tellus, when the or- 
bits are anti-aligned in the mean, and this result supports a similar 
conclusion following simulations illustrated in Figs.|4j|7] The sec- 
ular Vi resonance manifests itself as an upside-down V-like struc- 
ture (marked with yellow color) in the right-hand panel of Fig. [9] 
It is located between the 5c:2b and 3c: lb MMRs. Still, even for 
AG5b,c = 0, there is a large region above the 3c: lb MMR line, which 
is accessible for the HZ-stable orbits of Tellus. These results are in 
accord with shaded areas illustrated in Fig. 1 of |Malhotra & Minton] 



( |2008| ). Here, we found that mode II configurations (characterized 
by librations of AG5b,c around n for the neighboring trajectories), 
permit stable orbits of Tellus for almost entire observationally de- 
termined range of semi-major axes and eccentricities implying dy- 
namically stable configurations of the jovian sub-system. 



5 CONCLUSIONS 

The detection levels of extrasolar planets reach masses of super- 
Earths (i.e., Neptune-mass objects). The results of recent RV and 
transit surveys [Bouchy et al] ( [2008] ); |Bakos et al.| ( [2009] ), Udry 
2008 (an invited talk in Torun conference Extrasolar Planets in 
Multi-body Systems)] suggest that such small planets are common. 
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Figure 7. The left-hand panel. Color coded ^hz in the 5-plane for a\y = 2.3 au, ac = 4.6 au (model /o,7r). Gray crosses indicate unstable configurations of 
Jupiters in terms of (F) > 10. Nominal elements for A03b,c =0,71 are labeled with /q and 7;^, respectively. Dotted areas are for initial conditions implying 
librations of Acnb,c around (the right-hand half -plane of S) and around tt (the left-hand half-plane of S). The right-hand panel. The phase diagram calculated 
for a constant level of the total angular momentum (marked with gray curve in the 5-plane). 




Figure 8. Color coded coefficient of habitability, ^hz, illustrated in the 5- 
plane. Jovian semi-major axes are = 2.5 au, and ac =4.1 au, AC(3b,c = 
(model Vq, planets involved in the 2c: lb MMR). Gray crosses indicate 
strongly unstable configurations of the jovian sub-system, with (F) > 10. 
The tested configuration (model Vq, see Table[T} is marked with a circle. 
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Figure 9. Color coded coefficient of habitability, ^hz, illustrated in the 
(«b,«c) -plane. The left-hand panel is for A03b,c = tt, and the right-hand 
panel is for Acnb,c = 0. Eccentricities of the jovian planets are e\j = 
0.03, = 0.05 (models /o,7r). Gray crosses indicate chaotic configurations 
related to MMRs (labeled) of the b-c pair, (Y) > 10. 



Hence, we expect a growing interest in dynamical studies of ex- 
trasolar systems with terrestrial planets. Although we focus on the 
best analog of the Solar system ( [Gaudi et aLpOOSj ), in fact, our 
work concerns a whole class of multi-planet systems involving a 
hypothetical Earth-like object in the HZ of a low-mass star host- 
ing also large, jovian planets. To investigate the structure of the 
HZ, we adopted a simple quasi- analytic theory applicable to hier- 
archical configurations of the terrestrial planet and jovian compan- 
ions. Then the averaging principle can be applied, and the four- 
body problem may be further simplified in terms of the restricted 
model. Our quasi-analytic approach makes it possible to investi- 
gate the secular motion in an uniform way. Both the resonant and 
non-resonant orbits of primary bodies (Jupiter-mass planets) may 
be averaged out numerically. The secular equations of motion de- 
rived from the analytic theory, make it possible to reduce the CPU 
time by a few orders of magnitude. Hence, the method is useful 



to review the global features of the secular dynamics. In accord 
with the assumptions of the averaging principle, the quasi- analytic 
model can be also easily adopted to study long-term evolution of 
similar, dynamically scaled configurations, for instance the motion 
of jovian planets in multiple systems with brown dwarfs or sub- 
stellar companions. 

Regarding the OGLE-2006-BLG-109L system, we found that 
the structure of its HZ strongly depends on orbits of Jupiters, and 
not only on their semi-major axes and eccentricities, but also on 
the orbital phase and a libration mode of apsidal angle AG5b,c- From 
the recent theory we know ( [Michtchenko & Malhotra|2004| ), that 
apsidal librations around or 7i are generic orbital states of non- 
resonant, two-planet systems. In general, we found that the anti- 
aligned mode of the jovian orbits favors much larger fraction of 
stable orbits remaining entirely in the annulus of HZ between or- 
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bital radii of 0.25 an and 0.36 au, than are permitted by the aUgned 
configurations. Having in mind unconstrained parameters of the 
OGLE-2006-BLG-109L system, many scenarios are possible, even 
assuming that only one terrestrial planet resides in the HZ. The dis- 
astrous Vi secular resonance can be moved out of the HZ, through 
perturbations of additional smaller planets ( [Malhotra & Minton] 
|2008 ). They have shown that if such planets exist in the system then 
the secular dynamics of Earth-like planets would be very complex. 
Here, we also demonstrate that the dynamics are rich in the realm 
of the restricted four-body problem. 

Still, there remains an open question whether the creation of 
Earth-like bodies is possible in the OGLE-2006-BLG-109L sys- 
tem. We did preliminary simulations (Musielinski et al., in prepa- 
ration) with the help of MERCURY code ( |Chambers|1999| ), setting 
orbital parameters of the Jovian system consistent with the error 
bounds ( Gaudi et al. 20 08 ). We applied the common model of elas- 
tic coagulation of small protoplanets and planetesimals (see, e.g., 
|Raymond|2008| ). The initial distribution of planetary embryos in- 
volved ^ 50 Moon- sized objects distributed between 0.2 au and 
2 au. Curiously, as a typical outcome from 50 runs of 50 Myrs each, 
we obtained single 0.3-0.8 Earth-mass planet around 0.3 au and 
(sometimes) a second object of a similar mass beyond 0.6-0.8 au. 
The eccentricity distribution in the small sample of runs is basi- 
cally uniform, with a number of quasi-circular orbits. These results 
are encouraging but they must be confirmed by new, intensive sim- 
ulations. Nevertheless, we already found some evidence, that the 
terrestrial planets may emerge in the HZ of the OGLE-2006-BLG- 
109L system, and we have got a new motivation to investigate the 
effects of the Vi resonance in this region. 

In this work, we applied a few different techniques to study the 
dynamics of planetary systems. In particular, we apply the fast indi- 
cator (the symplectic MEGNO algorithm) to map the phase space 
of the system and to detect chaotic (unstable) configurations. The 
numerical integrations of the equations of motion are necessary to 
calibrate the integration time of the fast indicator and to test the 
quality of analytic expansions with the hep of semi-analytical av- 
eraging. To study the structure of the HZ, it is necessary to analyse 
simultaneously the short-term and the long-term behavior of the 
system. As we think, our study benefits from the synthesis of ana- 
lytical and numerical techniques. 
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